Multigrid solver for axisymmetrical 2D fluid equations 
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We have developed an efficient algorithm for steady axisymmetrical 2D fluid equations. The 
algorithm employs multigrid method as well as standard implicit discretization schemes for 
systems of partial differential equations. Linearity of the multigrid method with respect to 
the number of grid points allowed us to use 256 x 256 grid, where we could achieve solutions 
in several minutes. Time limitations due to nonlinearity of the system are partially avoided 
by using multi level grids(the initial solution on 256 x 256 grid was extrapolated steady 
solution from 128 x 128 grid which allowed using "long" integration time steps). The fluid 
solver may be used as the basis for hybrid codes for DC discharges. 



1. Introduction 

Further understanding of basic processes in gas dis- 
charges and non-equilibrium plasmas relies on com- 
parisons of experimental results and predictions of 
theoretical results which almost always have to be 
numerical calculations including self consistent cal- 
culation of the spatial profile of electric field. In 
all calculations it is critical to calculate the prop- 
erties of electrons which have large mobility and 
consequently ability to gain energy from the elec- 
tric field. Consequently electrons play a critical 
role in sustaining the plasma by gas phase ioniza- 
tion. At the same time electrons have non-local 
or non-hydrodynamic kinetics as their properties 
in the rapidly changing fields may not be defined 
uniquely by the local electric field. Thus, hybrid 
models were developed to take into account high 
energy electrons originating from the electrodes or 
created in very high electric fields close to the elec- 
trodes by using a Monte Carlo technique while the 
bulk of electros at low energies is accounted for by 
a fluid model. 

As the need to describe accurately more and 
more complex geometries transition from ID to 2D 
ad 3D systems becomes increasingly demanding in 
terms of computer time and complexity of equa- 
tions that may lead to numerical problems such as 
numerical diffusion and others. In particular as the 
grid becomes denser the computational demands in- 
crease as a square or cube of the number of grid 
points. Thus, in order to model realistic structures 
and complex geometries with proper relaxation of 
high energy electrons in space and time, one needs 
to develop special numerical procedures to handle 
the complex task. In this paper we describe one 
implementation of the multigrid technique which is 
far proved to be the leading contender for optimal 



treatment of 3D systems with a large number of grid 
points. 

2. Fluid equations and numerical algorithm 

Continuity equations for the electrons and the ions 
are [1] 



^ + V-(n e v e ) = S e , 

dn. 



+ V-(n p v p ) = S p . 



(1) 
(2) 



Momentum balance equations are 



cf> e = n e v e = -n e n e ~E- D e S/(n e ), (3) 
4> P = npV p = n p fi p E - D p V(n p ). (4) 

Electric potential is governed by the Poisson 
equation 

V 2 P = --(n p -n e ), (5) 
eo 

while the electric field is negative gradient of the 
potential: 

E = - VP. (6) 

In previous equations subscript index e (p) refers to 
the electrons (ions). Input parameters for the fluid 
equations are source terms S e and S p , mobilities 
Me(E) and /x p (E), diffusion coefficients D e (E) and 
D p(E) and are supposed to be known( e.g. in hybrid 
models they are provided from the lookup tables 
and Monte Carlo code). 

We will solve a set of equations ([H)-© for the 
azimuthal symmetry (/(r, z, G) = f(r, z)) and for a 
given set of boundary conditions. 

First we discretize system (P)-©: we split the 
domain (r, z) € [0, R] x [0, d] into rectangles with 
equidistant radial grid points ro, r±, . . . , rjsr r and non- 
equidistant axial grid points Zq,zi,..., zn z ■ We also 
have midpoints between the grid points r i+1 / 2 and 



z j+l/2- 



i+l/2 

In the following we use indices i,j,k for 



r-coordinate, z-coordinate and for the time, respec- 
tively. In order to allow long integration times the 
system must be discretized implicitly in time [2]. 
Continuity equations (pQ) and ([2]) are discretized "ba- 
ckward in time" 

nH 1 - n k ■ 

l ' J At ~ + (V.^)H 1 + (V^)H 1 = S hr (7) 

Momentum balance equations are of convection- 
diffusion type and it is convenient to discretize them 
by the Scharfetter-Gummel scheme [2]: 

<t>r,i+l/2,j = r ' 1 ^ 2 ' 3 ( n i,j G ( a l+l/2,j) ( 8 ) 



»z,i,j+l/2 



Azi 



■(«ijG(of J+1/2 ) (9) 
-n iii+ iG(a^. +1/2 )), 



where 



a 



Mr,i+l/2j 

a z,i,j + l/2 
ij+1/2 - _s 7) 



(Fi+ij - ^j), (10) 

r,,). (ii) 

z,i,j+l/2 

(with s = — 1 for the electrons and s = 1 for the 
ions) and 



F(x) 



G{x) 



x exp [x] 



exp [xj — 1 ' v ' exp [x] — 1 
Discretized form of the Poisson equation is 



+ 



8zj_i/ 2 Azj-iAzj 
(Ar) 2 fc+1 (Ar) 2 fc+1 



fc+i 



+ l 1 -2i) V '<-M + ( 1+ « 



V? 



fc+1 



2i 7 i+1 > j 



e /a \2 fc+1 i e / a \2 fc+1 n 

— n (Ar) n ei T -| n (Ar) n • , = 0, 

eo ' ,J eo 



(i = l,2,...,iV r -l;j = l,2,...,^ 2 -l), (12) 
while along direction r = is 

. (Ar) 2 (A^_ 1 + A^U fc+i , , T/ fc+i 



i 1/2AC, :A:, ) ^ + 



+ 



(AO 2 F fc+i (Ar) 2 fc+1 



z j-l/2^^-_i 



e / a \2 fc+1 1 e /a \2 fc+1 n 

— n (Ar) n ■ + — n (Ar) n ■ = 0, 
eo J eo 

(j = l,2,...,N z -l). (13) 



Discretized set of coupled equations connects 
concentrations and potential at time {k + l)At with 
values of concentrations and potential at time kAt. 
Since system of algebraic equations (f7|)- ([T3j) is non- 
linear we will solve it by applying the Newton- 
Raphson algorithm [3]. Linearizing the system by 
introducing a vector u = (n e , n p ,V) T and a correc- 
tion <5u of the same vector (during Newton- Raphson 
iterations), we obtain a set of linear equations for 
the correction vector: 

aij5ui-ij + bijSuij + CijSiii+ij (14) 
+dij5uij^i + eijSuij+i = fjj, 

(i = 0,l,...Ay-l;j = 1,2,...A+-1). 

where a denotes that a is a 3 x 3 matrix. Com- 
ponents of the matrices from the previous equation 
are long expressions and can be obtained straight- 
forwardly. Equation (|14p is computationally very 
demanding and one has to solve it very efficiently 
in order to get solutions of the fluid equations in a 
reasonable time. Having in mind that fluid solver 
is usually combined with Monte Carlo simulation in 
hybrid models, and that after every cycle of Monte 
Carlo simulation one solves the set of fluid equa- 
tions, the importance of the efficiency is evident. 

From known values of n e ,n p ,V at kAt we cal- 
culate the values at (A; + I) At by solving (fl~4"j) . Af- 
ter that we iterate (f7|)- ()13p through time until the 
steady concentrations are obtained. 

Since neither standard algorithms (e.g. Gauss 
elimination, LU decomposition) nor iterative algo- 
rithms (Gauss-Seidel, Jacobi, SOR)[3] are efficient 
enough, we have applied the multigrid method [4] 
to solve (fill) . The latter method takes 0(N) oper- 
ation for a linear system with N unknowns, while 
the former take at least 0(N 2 ). Typically 80 x 80 
or even coarser grids were used so far while we can 
use much finer grids. 

Moreover we have developed an improvement of 
the standard multigrid algorithm. The main idea 
consist of solving system (I14p on smaller discretiza- 
tion grid N r x N z , and then for the grid 2N r x 2N Z 
as an initial solution to take extrapolated station- 
ary solution from N r x N z and long integration time 
step. In such way we go to finer and finer grids, 
reaching the wanted number of grid points in the 
end (the finest one). 

The solution on N r x N z was obtained taking 
arbitrary initial n e ,n p and V and the initial time 
step was usually At\ = Ins (or any time for which 
the system (|14p is convergent). After each time 
integration step Ai&, new time integration step is 



Atk+i = nAtk (usually n = 5). It turns out that 
we can take greater and greater integration time 
steps as we approach the stationary solution. How- 
ever if tk + \ is too large (the norm of a correction 
vector in Newton-Raphson iteration exceeds some 
upper limit and the system is not robust for that 
time step), we go back and choose Atk+i = Ai&. In 
such way we approach the stationary solution very 
fast. We call our method multilevel prolongation 
method. 

3. Results 

On the basis of glow discharge measurements in ar- 
gon in cylindrical geometry with / = 920 fiA, V = 
255.9 V; d = 1.1 cm, R = 2.7cm, pd = 75 Pa ■ cm 
(p = 0.514 iorr)[6], we have have obtained an an- 
alytic form for the source terms S p (r, z) = 0.35 * 

|22 ^ (f) 2 exp(-50( f ) 8 ) 
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Changing the parameters a and b one may ob- 
tain different radial profiles which correspond to the 
cases of constricted and non-constricted discharge. 
Other input parameters are fx e = ^ m 2 V~ 1 s~ 1 

where p is the pressure given in torr; D e = lij^r 
with kT e = 1 eV; fi p was taken from [7]; D e = [i-p—f- 
0.026 eV. The boundary conditions are 



with kT p 

V(r,0) = 0; V(r,d) = 255.9 V; n e (r,0) 



n, 



(r,d) 



0; n p (r, 0) = n p {r,d) = 0; V(R, I) along the walls is 



interpolated linearly; n e 



n r 



along the walls. 



We have carried out calculations on a personal 
computer with processor Athlon 3200+ and 512 Mb 
of RAM memory. 

Number of Newton-Raphson iterations during 
the solving of nonlinear system was less than 10. 
Newton-Raphson iterations are stopped when the 
sum of the corrections was less than 10~ 6 . Time 
evolution are stopped when the relative change in 
densities become less than 10 -6 . The overall inte- 
gration time for the stationary solution was of the 
order of Is which is unimportant for our algorithm 
which advances through time by a factor n. 

Execution times are given in FigfTJ We can eas- 
ily see that the slope of the multigrid solvers (with- 
out (M) and with our acceleration (M+L)) is dif- 
ferent that for the iterative method which confirms 
different efficiencies of these methods. 

Figj2]-FigH] and Figj5j-Figl7] show the solution 
for a = 30,6 = 24 and a = 10, b = 1, respectively. 
We can see in FigJUFigJT] that the electron and ion 
densities are highly sensitive to the grid size. In the 
first case both densities increase with the grid, while 
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Figure 1: Execution times in seconds for different 
integration strategies. R denotes results obtained 
by iterative algorithm, M by multigrid, while L de- 
notes multilevel prolongation method. Grids have 



N r = N z . 




Figure 2: Source terms and densities for a = 30, 
b = 24 for grid size 64 x 64. 

in the second case the electron density decreases. 
The ion current to the cathode in the first case is 
/ = 0. 59771^4 and in the second case I = 0.035mA 
and are only very weakly dependant of the grid size. 

4. Conclusions 

We have developed an efficient method for 2D fluid 
equations in cylindrical geometry. We have applied 
the multigrid to our knowledge for the first time 
to solve a system of equations for an obstructed dc 
glow discharge and with implementation for multi 
level prolongation approach. Tremendous accelera- 
tion in respect to the iterative methods is obtained 
by applying multigrid method combined with mul- 
tilevel prolongation method. From the solution we 
may conclude importance of the grid size: is not so 
important for the ion flux to the cathode, but for 
densities it is important. 

A question which arises is how fine grid may be 
imposed by the physics of the problem. The an- 
swer should be probably not finer than the mean 



Figure 3: Potential and electric field for a = 30, 
b = 24 for grid size 64 x 64. 
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Figure 4: Potential and electric field for a = 30 
b = 24. 

free path for the ionization which may vary with 
the external parameters of the discharge. 
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Figure 5: Source terms and densities for a = 10, 
b = 1 for grid size 64 x 64. 




Figure 6: Potential and electric field for a = 10, 
6 = 1 for grid size 64 x 64. 
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Figure 7: Potential and electric field for a = 10 
6=1. 



